# -*- coding: utf-8 -*-
"""
Created on Thu Jul 28 08:18:36 2022

@author: Derek Rodriguez Pacheco, PhD Candidate
IUSS Pavia
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import subprocess
from scipy.optimize import minimize


betaf= 1
pinchX= 0.5
pinchY= 0.01

params = [betaf,pinchX,pinchY]



def HystModelFit(par):

    k1= 0.5
    k2= 0.1
    sigAct= 3.4
    epsSlip= 0
    epsBear= 0
    rBear= 1
    
    s1p= 2.5
    e1p= 6
    s2p= 4
    e2p= 25
    s3p= 15
    e3p= 20
    
    damage1= 0
    damage2= 0
    beta= 0
    
    s1n= s1p*45
    e1n= np.copy(e1p)
    s2n= s2p*45
    e2n= np.copy(e2p)
    s3n= s3p*45
    e3n= np.copy(e3p)
    
    m= 970                           ## Mass in N*s2/m or kg
    g= 9.81                          ## Acceleration of gravity in m/s2
    
    adjfactors= np.array(pd.read_csv('Adjustment_Factors_Displacement_CM.txt', header= None, delim_whitespace= True)[0])
    indexesb= np.array(pd.read_csv('Indexes_Beginning.txt', header= None, delim_whitespace= True)[0])
    indexesf= np.array(pd.read_csv('Indexes_End.txt', header= None, delim_whitespace= True)[0])
    
    ## Model's parameters
    
    
    
    width= 800
    height= 1076
    
    testnumber= [13,14,15,16,17,18,19,20,24,25,26,27,28,40,41]
    
    ##
    
    datacm= {}
    # acceltcm01= {}
    # acceltcm02= {}
    accelcm= {}
    timecm= {}
    timecmc= {}
    
    data= {}
    datad= {}
    accelda= {}
    # accelda1= {}
    # accelro= {}
    accelbase= {}
    dispda= {}
    dispdacm= {}
    timearray= {}
    energyarray= {}
    energyarraycm= {}
    timetotal= []
    acceltotal= []
    dispnet= {}
    dispnetcm= {}
    # timestart= []
    
    for i in range(15):
        datacm['%d' % (i+1)]= np.array(pd.read_csv('Data_Center_Mass\\accelAndDispl_15Hz\\%d_AccelAndDispl_15Hz.txt' % testnumber[i], header= None, delim_whitespace= True))
        accelcm['%d' % (i+1)]= np.array(datacm['%d' % (i+1)][:,4])[indexesb[i]:indexesf[i]+1]
        timecm['%d' % (i+1)]= np.array(datacm['%d' % (i+1)][:,0])[indexesb[i]:indexesf[i]+1]
        timecmc['%d' % (i+1)]= timecm['%d' % (i+1)]-timecm['%d' % (i+1)][0]
    
   
    data['1']= np.array(pd.read_csv('13_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['2']= np.array(pd.read_csv('14_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['3']= np.array(pd.read_csv('15_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['4']= np.array(pd.read_csv('16_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['5']= np.array(pd.read_csv('17_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['6']= np.array(pd.read_csv('18_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['7']= np.array(pd.read_csv('19_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['8']= np.array(pd.read_csv('20_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['9']= np.array(pd.read_csv('24_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['10']= np.array(pd.read_csv('25_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['11']= np.array(pd.read_csv('26_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['12']= np.array(pd.read_csv('27_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['13']= np.array(pd.read_csv('28_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['14']= np.array(pd.read_csv('40_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    data['15']= np.array(pd.read_csv('41_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
    
    datad['1']= np.array(pd.read_csv('13_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['2']= np.array(pd.read_csv('14_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['3']= np.array(pd.read_csv('15_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['4']= np.array(pd.read_csv('16_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['5']= np.array(pd.read_csv('17_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['6']= np.array(pd.read_csv('18_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['7']= np.array(pd.read_csv('19_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['8']= np.array(pd.read_csv('20_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['9']= np.array(pd.read_csv('24_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['10']= np.array(pd.read_csv('25_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['11']= np.array(pd.read_csv('26_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['12']= np.array(pd.read_csv('27_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['13']= np.array(pd.read_csv('28_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['14']= np.array(pd.read_csv('40_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    datad['15']= np.array(pd.read_csv('41_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
    
    for i in range(15):
        accelda['%d' % (i+1)]= np.array(data['%d' % (i+1)][:,4])
        accelbase['%d' % (i+1)]= np.array(data['%d' % (i+1)][:,1])
        dispda['%d' % (i+1)]= np.array(datad['%d' % (i+1)][:,3])
        dispnet['%d' % (i+1)]= np.array(datad['%d' % (i+1)][:,4])
        dispdacm['%d' % (i+1)]= np.array(datad['%d' % (i+1)][:,3])*adjfactors[i]
        dispnetcm['%d' % (i+1)]= np.array(datad['%d' % (i+1)][:,4])*adjfactors[i]
        timearray['%d' % (i+1)]= data['%d' % (i+1)][:,0]-data['%d' % (i+1)][0,0]
        energyarray['%d' % (i+1)]= [np.trapz(accelda['%d' % (i+1)][:j]*m,dispda['%d' % (i+1)][:j]/1000) for j in range(len(accelda['%d' % (i+1)]))]
        energyarraycm['%d' % (i+1)]= [np.trapz(accelcm['%d' % (i+1)][:j]*m,dispdacm['%d' % (i+1)][:j]/1000) for j in range(len(accelcm['%d' % (i+1)]))]
        
    cc= 0
    
    for i in range(15):
        for j in range(len(timearray['%d' % (i+1)])):
            if i==0:
                timetotal.append(timearray['%d' % (i+1)][j])
            else:
                timetotal.append(cc+timearray['%d' % (i+1)][j])
            acceltotal.append(accelbase['%d' % (i+1)][j])
        cc= timetotal[-1]
        if i!=14:
            for j in range(int(round(15/timearray['%d' % (i+1)][1]))):
                if i==0:
                    timetotal.append(timearray['%d' % (i+1)][-1]+(j+1)*timearray['%d' % (i+1)][1])
                else:
                    timetotal.append(cc+(j+1)*timearray['%d' % (i+1)][1])
                acceltotal.append(0)
            cc= timetotal[-1]
    
    #### SDOF Analysis
    
    # for i in range(15):
    np.savetxt('TH/Time_Total.txt', timetotal)
    np.savetxt('TH/Acceleration_Total.txt', acceltotal)
    
    recordsa= 'TH/Acceleration_Total.txt'
    recordst= 'TH/Time_Total.txt'
    
    dts= 0.00195
    mass= m/1e6
    
    dur= timetotal[-1]
    
    fh= open('FMdata.tcl', 'w')
    
    fh.write(
        'set FMfilet "%s";\nset FMfilea "%s";\nset dtg %f;\nset durs %f;'
        '\nset M %f;'
        '\nset s1p %f;\nset e1p %f;\nset s2p %f;\nset e2p %f;\nset s3p %f;'
        '\nset e3p %f;\nset pinchX %f;\nset pinchY %f;\nset damage1 %f;'
        '\nset damage2 %f;\nset beta %f;' 
        '\nset s1n %f;\nset e1n %f;\nset s2n %f;\nset e2n %f;\nset s3n %f;\nset e3n %f;'
        '\nset width %d;\nset height %d;'
        '\nset k1 %f;\nset k2 %f;\nset sigAct %f;\nset betaf %f;\nset epsSlip %f;\nset epsBear %f;\nset rBear %f;' % (recordst, recordsa, dts, dur, mass, s1p, e1p, s2p, e2p, s3p, e3p, par[1], par[2], damage1, damage2, beta, s1n, e1n, s2n, e2n, s3n, e3n, width, height, k1, k2, sigAct, par[0], epsSlip, epsBear, rBear))
        # % (recordst, recordsa, dts, dur, mass, s1p, e1p, s2p, e2p, s3p, e3p, pinchX, pinchY, damage1, damage2, beta, s1n, e1n, s2n, e2n, s3n, e3n, width, height))
    fh.close()
    subprocess.call('OpenSees MDOF_CalibrationMH_Grav.tcl', shell=True)
    
    ####
    
    sdofd= {}
    sdofa= {}
    sdoft= {}
    sdoff= {}
    sdofe= {}
    
    indexes= [0,22058,44116,66183,88163,110144,132124,154103,176104,198072,219183,240296,261597,283035,304762,319058]
    timei= [0,43.014,86.028,129.057,171.918,214.781,257.642,300.501,343.403,386.241,427.407,468.579,510.116,551.919,594.286]
    
    for i in range(15):
        sdofd['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/DispT.out', header= None, delim_whitespace= True)[1][indexes[i]:indexes[i+1]])*0.001
        sdofa['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/AccelT.out', header= None, delim_whitespace= True)[1][indexes[i]:indexes[i+1]])
        sdoft['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/AccelT.out', header= None, delim_whitespace= True)[0][indexes[i]:indexes[i+1]])-timei[i]
        sdoff['%d' % (i+1)]= sdofa['%d' % (i+1)]*0.001*m*-1
        sdofe['%d' % (i+1)]= [np.trapz(sdoff['%d' % (i+1)][:j],sdofd['%d' % (i+1)][:j]) for j in range(len(sdoff['%d' % (i+1)]))]
        
    rtzero= [14365,14365,14373,14287,14287,14287,14287,14307,14276,13933,13420,13610,13743,14035,14296]
    
    ErrorEn = 0.0
    ErrorDisp = 0.0
    ErrorFor = 0.0
    
    maxE = np.zeros(15)
    maxD = np.zeros(15)
    maxF = np.zeros(15)
    
    for i in range(1,16):
       maxE[i-1] = max(np.absolute(energyarraycm['%d' % i]))
       maxD[i-1] = max(np.absolute(dispdacm['%d' % i]))
       maxF[i-1] = max(np.absolute(accelcm['%d' % i]))
       
    #Enmax = max(maxE)
    #Dispmax = max(maxD)
    #Formax = max(maxF)
            
            
    for i in range(1,16):
        Net = len(energyarraycm['%d' % i])
        Nen = len(sdofe['%d' % i][:rtzero[i-1]])
        Ne = min([Net,Nen])
        for j in range(Ne):
            ErrorEn = (energyarraycm['%d' % i][j]/maxE[i-1]-sdofe['%d' % i][j]/maxE[i-1])**2
            ErrorDisp = (dispdacm['%d' % i][j]/maxD[i-1]-sdofd['%d' % i][j]*1000/maxD[i-1])**2
            ErrorFor = (accelcm['%d' % i][j]*m/maxF[i-1]-sdoff['%d' % i][j]/maxF[i-1])**2
            
    wDisp = 0.5
    wFor = 0.1
    wEn = 0.4
    
    TotError = wDisp*ErrorDisp+wFor*ErrorFor+wEn*ErrorEn
    
    return TotError

HystModelFit(params)
    
res = minimize(HystModelFit,params,method = 'Nelder-Mead',options={'disp': True})  

res.x


